States without linear counterpart in Bose-Einstein condensates 
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We show the existence of stationary solutions of a 1-D Gross-Pitaevskii equation in presence of 
a multi-well external potential that do not reduce to any of the eigenfunctions of the associated 
Schrodinger problem. These solutions, which in the limit of strong nonlinearity have the form of 
chains of dark or bright solitons located near the extrema of the potential, represent macroscopically 
excited states of a Bose-Einstein condensate and are in principle experimentally observable. 

PACS numbers: 03.65.Ge, 03.75.Fi, 47.20.Ky 



I. INTRODUCTION 

Bose-Einstein condensation (BEC) of weakly interact- 
ing atomic gases [Q strongly motivates the study of the 
Gross-Pitaevskii equation (GPE), 



h 2 , 

2m 



U Q \^(x,t)\ 2 + V(x) 



*(sc,t) = 0, 



(1) 

a mean-field Schrodinger equation with local cubic non- 
linearity. Of particular interest are the non ground-state 
stationary solutions of the GPE ||] which represent 
macroscopically excited states of the condensate. Vor- 
tices have been recently observed in two- Q| or one- 
component H condensates and are also invoked as a su- 
perfluidity breaking mechanism |^|. Phase engineering 
optical techniques have allowed to generate dark solitons 
in atomic gases with positive scattering length |Q ||. 

Vortices and solitons observed in recent experiments 
are examples of excited states with linear counterpart, i.e. 
stationary solutions of the GPE which can be obtained 
as a deformation of eigenstates of the corresponding lin- 
ear Schrodinger equation [u|. However, the GPE may 
also admit stationary solutions without linear counter- 
part. In a discretized version of the GPE, also known 
as discrete self trapping equation, the existence and sta- 
bility of solutions without linear counterpart has been 
studied at various discretization orders pd[ |. In particu- 
lar, the appearance of self trapping stationary states in 
the dimer case, which mimics a double-well system, has 
been widely investigated in connection with the evolution 
of wave packets [O, 13, 14]. Recently, a set of stationary 
solutions without linear counterpart has been discovered 
also in the continuous case, namely the exactly solvable 
1-D GPE with periodic boundary conditions and zero ex- 
ternal potential Jl(J. These states break the rotational 
invariance of the associated linear problem. 
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In this paper, we show the existence of stationary so- 
lutions without linear counterpart of a 1-D GPE in pres- 
ence of a multi-well external potential. In the limit of 
strong non linearity, these solutions assume the form of 
chains of dark or bright solitons located near the extrema 
of the potential and in general break the symmetry of the 
external potential. 

Our analysis is of direct interest for BEC experiments 
where atomic gases can be confined in arbitrarily tailored 
magnetic or optic traps. As a case study, we investigate a 
GPE representing a quasi 1-D Bose-Einstein condensate 
confined in a double- well trap described by the potential 



V(x) = m 2 7 4 x 4 
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In Section III we describe all the zero-, one-, and two- 



soliton solutions of this model in an analytical way valid 
in the limit of strong non linearity. In Section ^ by 
means of numerical simulations we find the exact shape 
of these states and study their evolution in the linear 
limit reached when the number of particles in the conden- 
sate, N, vanishes. We consider both the cases of conden- 
sates with positive or negative scattering length. Shape 
and energy of the corresponding stationary solutions are 
shown in Fig.s 1, 3 and 2, 4, respectively. Their stability 
properties are discussed in Section [v|. 



II. STATIONARY SOLUTIONS 

Let us first review some general properties of the sta- 
tionary solutions of the GPE that reduce, in the limit of 
vanishing nonlinearity, to the eigenfunctions of the asso- 
ciated Schrodinger equation 



2m 



V 2 + V{x) - £ r , 



4> n {x) = 0, n = 0, 1, . 



(3) 
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In U we have shown that for any finite value of the chem- 
ical potential /x there exists a set of stationary solutions 
of the GPE, Of fj, n (x , t) = exp (— i/zt) ^^(x), which have 
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FIG. 1: Zero-, one-, and two-soliton stationary solutions of the repulsive GPE with the symmetric double-well potential (|2|) for 
different values of the normalization N. For comparison, the functions are shown scaled by y/N. The vertical solid and dashed 
lines indicate the double-well maximum and minima, respectively. The degenerate states obtained by changing ip{x) — > ip(—x) 
are not reported. The results have been obtained with the following parameters: m — 3.818 x 10 -26 Kg, lo = 12.75 Hz, 
7 = 10 9 Kg-im-is-i, U = 1.1087 x 10" 41 Jm. 
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FIG. 2: As in Fig. |l| for the one- and two-soliton solutions in the attractive case Uq = —1.1087 x 10 41 Jm. 



limit ipf.nix) \\ipnn\\ 1 ~ > <l>n(x) when n — > £ n . The pa- 
rameter \i ranges in the interval [£ n , +oo) for Uq > and 
in the interval (— oo,£ n ] for Uq < 0. In both cases, the 
number of particles in the state i/Vm N n ([i) = ||VVn|| 2 , 
vanishes for fj, — > £ n . In other words, the linear limit is 
reached for a vanishing number of particles in the con- 
densate. 

In the 1-D case, asymptotically exact expressions for 
the GPE stationary solutions with linear counterpart are 
known also in the opposite limit of strong nonlinearity. 
For — > ±oo, depending on the sign of Uq, these solu- 
tions assume the form of chains of dark or bright solitons 
P). More specifically, in the repulsive case Uq > the 



solution with n — nodes assumes the zero-soliton shape 



while for n > 1 nodes we obtain asymptotic solutions 
with n dark solitons 



^^n{x) -> ty^x) Y[ tanh ( - x k ) ) . (5) 

fc=i ^ 



In the attractive case Uq < 0, for /i — > — oo the solutions 



3 



with n > nodes give rise to n + 1 bright solitons 




sech 



fc=0 



(x - x k ) 



(6) 



In the functions (|||j) with two or more solitons, the 
solitons do not overlap, i.e. the distance between their 
centers Xk is much larger than the dark-soliton width 
h/y/mfj, or the bright-soliton width H/ ^/—2m^i 



Note 

that any stationary solution is invariant under a global 
phase change and we do not consider this trivial degen- 
eracy. 

The stationary solutions of the GPE, for fi fixed, are 
the critical points of the grand-potential functional 



n[4>] = 



:W*)| S 



dec. 



(7) 



Since for \fi\ large the GPE solutions with linear counter- 
part assume the form (fi]||) with some specified centers 
{xfc}, we look for more general multi-soliton solutions 
in which the soliton centers may assume different val- 
ues. The allowed {xk} can be determined by substituting 
the expressions ((§](]) in (^) and extremizing the resulting 
function fl({xk})- 



III. ZERO-, ONE-, AND TWO-SOLITON 
SOLUTIONS IN A DOUBLE- WELL 

Zero-soliton solutions exist only in the repulsive case 
Uq > and are given by Eq. (Q). For [i sufficiently 
large, we have a node-less state which extends over the 
entire double well (column 1 of Fig. 1). If fj, is smaller 
than the barrier height uj^/Aj 4 , this state vanishes in 
the barrier region where V(x) > /i. In this case, since 
ip = is a trivial solution of the GPE, we could expect 
also two other stationary solutions of the form ip(x) — 
y{pi — V(x)) JUq in one of the two wells and ip(x) = 
elsewhere (column 2 of Fig. 1 and symmetric partner 
ip(—x)). The new solutions break the symmetry of V 
and must disappear in the linear limit. They correspond 
to the self-trapped states studied in JD], |l2|, [l3|, [l4| . 

One-soliton solutions are described by Eq. (pJ) with 
n = 1 in the repulsive case and Eq. (^) with n — in 
the attractive one. The corresponding grand-potential 
becomes a function of the soliton centers Xi or xq, re- 
spectively. For \n\ sufficiently large, the width of the 
solitons is very small and the dependence of the integral 
(f?J) on x\ or xq is due only to the term F|'0| 2 . The dark 
is constant with a hole in x\ so 
V(x\). The bright soliton density 



soliton density \tp^i 
that VL(x\) ~ const 

| "0^0 1 2 is different from zero only in proximity of xq and 
^I(xq) ~ const + V(xq). In both cases we have three one- 
soliton solutions corresponding to the three extrema of 



the external potential. The soliton may be found in the 
maximum (column 3 of Fig. 1 and column 1 of Fig. 2) 
or in one of the two minima (column 4 of Fig. 1 and col- 
umn 2 of Fig. 2 and symmetric partners ip(—x)) of the 
double-well. The two so lutions wi th the soliton centers 
in ±x m , where x m = ^/w 2 /2m7 4 , break the symmetry 
of V(x) and do not have linear counterpart. 

In the repulsive case, two-soliton solutions are de- 
scribed by Eq. (^) with n = 2 and the grand-potential 
becomes the two-variable function a^)- When the 

distance between the soliton centers is much larger than 
their width, we have Cl(xx,X2) — f2(xi) + f2(x2). In the 
region x\ < x%, Q has a maximum in (— x m , x m ) and two 
saddle points in (0, x m ) and (—x m ,0). We assume that 
x m h/y/mjl. The stationary solution corresponding to 
the maximum of f2 is shown in column 5 of Fig. 1. Those 
corresponding to the two saddle points (column 6 of Fig. 
1 and symmetric partner ip(—x)) break the symmetry of 
V and must disappear in the linear limit. 

Other extrema of f2 can be found when the centers of 
the two dark solitons are into the same well. In fact, 
when both x\ and X2 tend to x m , or —x m , the value 
of Q(xi,X2) ~ const — V(xx) — ^(^2) increases until 
\x\ ~ x 2 \ ^> H/ \Jfn4L. When the distance \x\ — x 2 \ be- 
comes comparable with the soliton width, the two den- 
sity holes in IVV2I 2 begin to merge and the norm of '0/i2 
increases. This implies that f2 decreases for \x\ — x%\ — > 
since, at least for fi sufficiently large, S! ~ — A* 1 1 1 1 2 • As 
a consequence, has two maxima in (x m — S, x m + S) 
and (— x m — (5, — x m + 5) with 28 > h/ \jm\x. The corre- 
sponding solutions (column 7 of Fig. 1 and symmetric 
partner x)), break the symmetry of V(x) and do not 
have linear counterpart. 

In the attractive case the situation is more compli- 
cated. The bright solitons in the stationary solutions 
with linear counterpart given by Eq. (|^) are multiplied 
by a phase factor which is alternatively +1 and — 1. In 
general, we can expect bright solitons with arbitrary rel- 
ative phases since each sech function is, for fi — > —00, 
solution of the GPE and this equation is invariant un- 
der a global phase change. Restricting to real solutions, 
in the two-soliton case we have to consider the following 
possibilities 



. 1 V-2m^ 
sech I (x — Xq) 



, , y/^2mp 
±sech I (x — xi) 



(8) 



The functions ^.^(xo, x\) obtained by inserting these ex- 
pressions in (Q) present, in analogy with the repulsive 
case, a minimum in (— x m ,x m ) and two saddle points in 
(0,x m ) and (— x m ,0). The stationary states correspond- 
ing to the minimum of f2 ± (xo, x\) are shown in columns 
3 and 5 of Fig. 2. Those corresponding to the two saddle 
points (columns 4 and 6 of Fig. 2 and symmetric part- 
ners V ,± ( — x )) break the symmetry of V and do not have 
linear counterpart. 



4 




io 2 io 3 io 4 io 5 

N 



FIG. 3: Single-particles energies E/N for the states shown 
in Fig. | as a function of N. The numbers correspond to 
the columns of Fig. The curve 8 corresponds to the anti- 
symmetric partner of 5. 
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FIG. 4: As in Fig. ^| for the states shown in Fig. |[ Curves 
8 and 9 correspond to the stationary solutions, not shown in 
Fig. |j| having as linear counterpart the Schrodinger eigen- 
functions with 2 and 3 nodes, respectively. 



IV. NUMERICAL SOLUTIONS WITH AN 
ARBITRARY NUMBER OF PARTICLES 

Now we compare the zero-, one-, and two-soliton solu- 
tions discussed above with the results of numerical sim- 
ulations. We use a numerical algorithm based on a stan- 
dard relaxation method for partial differential equations 
p"5[ . The success of this method is crucially based on the 
quality of the trial functions used to start the relaxation. 
For \n\ very large, good trial functions are represented 
by the multi-soliton functions with the soliton centers 
determined as above. The relaxed solutions can be then 
used as trial functions for a new simulation with a smaller 
value of |/^| . By changing /i sufficiently slowly, one can fol- 
low the evolution of the stationary states until they reach 
the linear limit, if it exists, or the point where they disap- 
pear. Figures [I] and |2| show in the repulsive and attractive 
cases, respectively, the states obtained in this way for dif- 
ferent values of their norm N . For any node index n, the 
linear limit is reached when N = N n (fi) — » 0. All the so- 
lutions that break the symmetry of the external potential 
disappear for N smaller than a critical value. However, 
there also exist solutions without linear counterpart that 
preserve this symmetry. An example is shown in the first 
column of Fig. || which corresponds to a bright soliton 
at the center of the barrier. 

In Figs. | and | we show the single particle energies 
for the same states of Figs, [l] and || as a function of 
N. From these figures it is evident the generation of 
solutions without linear counterpart as N is increased. 
In the case of the attractive GPE, the stationary solution 
which for N large is fully localized into one of the two 
wells (second column of Fig. ||) is, when it exists, the 
state of minimal energy. Therefore, the nature of the 
mean-field ground state changes as a function of N and 
this suggests the existence of a quantum phase transition 
in the corresponding exact many-body system. 

The generation of stationary states without linear 
counterpart can be understood in terms of bifurcations of 
superpositions of Schrodinger eigenstates. In the follow- 
ing we discuss an analytical example valid when the zero 
point energy of each isolated well, ^fi2w, is much smaller 
than the barrier height, w 4 /47 4 , i.e. for oj 3, ^> I. Let 
us consider stationary solutions of the GPE of the form 

ip(x) = VN [a Q xo(x - x m ) + b Q xo(x + x m )] , (9) 



On the other hand, due to the gradient term in (^) , we 
have a different behavior of f2 + and fl~ when both the 
soliton centers xo and x\ move toward the minimum of 
the same well. In fact, fi + does not present new extrema 
while Q~ has two minima in (x m — 5, x m + 5) and (— x m — 
S, —x m + 8) with 2(5 > hj ^/-2ro/i. The corresponding 
solutions (column 7 of Fig. 2 and symmetric partner 
(— x)), break the symmetry of V(x) and do not have 
linear counterpart. 



where Xn(x) are the eigenfunctions of the Schrodinger 
problem with harmonic potential ^m(2Lu) 2 x 2 and a§ + 
6q = 1. Since the state (Q) is normalized to TV, for it to be 
a stationary solution of the GPE we have to extremize the 
energy functional E[ip] = fiN . Up to exponentially 

small terms we get 

E(b ) ~ b Jl-b? ) + sign(Uo)^- {l + 2b 4 - 2bl) , 

(10) 



5 



where 



N 



LU 

'Jvf 



^h 3 u>/mUl (11) 



For N <C Nq, E(b ) has a minimum for 6 — 2 2 and a 
maximum for bo = — 2~2. These extrema correspond to 
the lowest energy symmetric and anti-symmetric linear 
states (columns 1 and 3 of Fig. [j] and columns 3 and 
5 of Fig. ||). If Uq > 0, for N ~ N the maximum at 
bo = — 2~2 bifurcates in a minimum and a maximum 
which, increasing N, moves to bo = 0. This describes 
the birth of the state in the second column of Fig. |l| 
and its subsequent localization in the right well (do = 1)- 
If Uq < a similar result is obtained with maxima and 
minima exchanged (see column 2 of Fig. ^|). Generation 
of other states can be obtained by considering superpo- 
sitions more complicated than rt9|). 



V. STABILITY OF STATIONARY SOLUTIONS 

In this Section we discuss the stability of the stationary 
states described above. We start with a linear stability 
analysis. Consider the linearization of Eq. ([[]) for a small 
change of its solution 



dt 



2m 



V 2 + V(x)+2U \V\ 



5* + C/ * 2 (5**. 

(12) 



We are interested to evaluate the evolution of the varia- 
tion 5^ of a stationary solution By writing 



* + = * (U „(x,t) + 5*(a;,t) 
= e-*" t [V7m(*) + <ty( 



according to Eq. (^2() the variation 
conjugated 5(jf are determined by 



lh dl 



where 



V 



Uofpfj 
-V 



fin 



(x,t)], (13) 
and its complex 



(14) 



2V = ~^V 2 + V(x) + 2U \^ n \ 2 - 
2m 

The solution of Eq. (fL4f) can be written as 



(x,t) 

( X ,ty 



k 



i\ k t/h ( fk(x) 
9k(x) 



(15) 



(16) 



where Afc and (f k ,g k ) are the eigenvalues and the eigen- 
vectors of the linearization operator 



^ fin 

-u r u 



-T>„ 



fk 
9k 



fk 

9k 



(17) 



and the coefficients c k are fixed by the initial condition 
<ty(x,0). Multiplying @ by {f* kl -g* k ) and integrating 
over space, we get 



J [f£Dfmfk + flfcf/in 



(IM 



da;. 



(18) 



The linearization eigenvalues are real and Eq. (|1J) ad- 
mits quasiperiodic solutions for an arbitrary -0 M „ Jl6| . 
This proves the linear stability of the stationary solu- 
tions. 

Lyapunov stability of all the states * in the neighbor- 
hood of a stationary solution p„ is a mathematically 
stronger concept of stability and certainly more relevant 
from an experimental point of view. This kind of stabil- 
ity was previously studied in the case of a lattice model 
which reduces to the GPE in the continuum limit |l7| . In 
that paper it was shown numerically that the maximum 
Lyapunov exponent associated to a discretized version of 
Eq. (12) vanishes when the initial state ^(x,0) is suf- 
ficiently close to one of the stationary states. A similar 
analysis can be pursued in the case of the GPE by sim- 
ulating a huge finite dimensional system, namely that 
obtained by applying a finite difference scheme to the 
partial differential equation ([!]). Of course, the large but 
finite number of degrees of freedom used in the simulation 
sets a limit to the maximum time at which the properties 
of the infinite dimensional system corresponding to the 
GPE are correctly represented. We will report on this 
elsewhere. Here we note that the gained scenario is con- 
sistent with Kuksin theory |Q which asserts that Eq. 
( |l2| ) admits iV-dimensional invariant tori, deformation 
of (pi]), in a finite neighborhood of any stationary state 
whose linearization spectrum satisfy non-degeneracy and 
non-resonance conditions. 

The spectrum of the linearization operator is useful 
also for discussing the stability of the stationary states 
under the effect of a dissipative perturbation. The grand- 
potential (^) evaluated for a state of the form ([l3|), up to 
the second order in the variation 8cj> gives 



where 



s 2 n 



(19) 



da; 



i> [V^* + Uo% n H\ dx. (20) 
By using Eqs. ( |l4| ) and (|l6| ) and the sum rule 



E 



Cke 



-iXkt/h 



f k = J24e tXkt/h gl 



we get 



<5 2 fi = i$> fe | 2 A fe (||/ fe f 



\9k\ 



(21) 



(22) 
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Therefore, a stationary solution -0 M „ is a local minimum 
of the grand-potential functional if and only if for any A; 
we have 
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x k (\\fk\\ 2 -\\g k \\ 2 )>o. 



(23) 



To verify the disequalities (|23|), we have solved numer- 
ically the eigenvalue problem (y_7|) by representing the 
linearization operator with a finite difference scheme. In 
the repulsive case Uq > 0, the condition (^3|) is fulfilled 
only by the state in the first column of Fig. 1. In the 
attractive case Uq < 0, no one of the states shown in Fig. 
2 satisfies ( p3| ) . This can be explained observing that the 
grand-potential evaluated at a stationary state is 



1 



U / |V7m| da;. 



(24) 



Thus, ii Uq < 0, fl assumes the minimal value for the 
trivial solution ip = 0. These results have a certain in- 
terest on the stability of a physical condensate in which 
a dissipative dynamic is introduced by the coupling with 
the environment degrees of freedom. Eventually the sys- 
tem will converge to a local minimum of Q. For attractive 
interaction, this implies the disappearance of the conden- 
sate. An estimate of the characteristic lifetimes has been 
given in the case of a vortex state (20) . 

We conclude our stability analysis by considering the 
short time behavior of the stationary states under the 
action of an initial finite deformation. A similar analy- 
sis has been considered in 0] to check the stability of 
the solutions found in [n0|. The authors of |21 stud- 



ied the evolution of stationary states initially perturbed 
with a stochastic noise. In our case, the stationary solu- 
tions corresponding to solitons located near the extrema 
of the double- well potential should have great sensitivity, 
specially in the case of unstable extrema, to symmetry 
breaking perturbations. Here, we consider the evolution 
of shifted stationary states, i.e. we solve Eq. ([[]) with 
the initial condition ^(x, 0) = i]j^ n {x — Ax). The numer- 
ical simulations have been performed with the improved 
Crank-Nicholson scheme introduced in JItJ which pro- 
vides an accurate conservation of the constants of mo- 
tion of (|l|), namely norm and energy. As an example, 
we describe the evolution of the states of column 4 of 
Fig. 1 and 2. Note that these are states without linear 
counterpart. 

The state with one dark soliton has a rather simple 
evolution. The qualitative shape of the state does not 
change but the soliton oscillates around the minimum 
x = x m of the right well. The position of the soliton 
center, ii, as a function of time is shown in the upper 
panel of Fig. || for Ax — 0.3 //m. The amplitude of 
the oscillations is very small in the case considered and 
increases by increasing Ax. For a shift Ax sufficiently 
large the soliton can jump between the two wells. 

The dynamics of the state with two bright solitons is 
more complicated, see lower panel of Fig. H. Initially the 
soliton at the center of the barrier moves toward that 




FIG. 5: Time evolution of the soliton centers for initially 
perturbed stationary states ^(x, 0) = ippm{x—Ax) with Aa; = 
0.3 fim. In the upper panel ip^n is one of the states in column 
4 of Fig. 1 while in the lower panel one of those in column 4 
of Fig. 2. 



located inside the right well, the latter being essentially 
at rest. When the distance between the two solitons be- 
comes comparable to their width the oscillations of the 
solitons inside the same well turn out to be correlated. 
The solitons do not cross each other. When the poten- 
tial energy of the barrier, T^|4'| 2 , is sufficiently reduced 
by the negative interaction energy, if/ol^l 4 , the soliton 
which was originally at the center of the barrier can jump 
into the left well. This inter-well dynamics is obtained 
also for values of the initial shift smaller than Aa; = 0.3 
fim which is the case shown in Fig. By decreasing 
Ax, the initial falling of the soliton at the center of the 
barrier into the right well (left well for Ax < 0) is slowed 
down. 

The results shown in Fig. |^ can be generalized to dif- 
ferent kinds of perturbations, e.g. stochastic noise, mod- 
ification of the parameters of the external potential. De- 
tails will be reported elsewhere. 



VI. CONCLUSIONS 

We have shown that in presence of an external poten- 
tial a 1-D GPE can admit stationary solutions without 
linear counterpart. Their existence is strictly connected 
to the multi-well nature of the potential. In the dou- 
ble well example discussed here, these solutions disap- 
pear in the limit ui — > when the potential assumes the 
shape of a single quartic well. For a piece-wise constant 
double-well, the stationary states here discussed analyt- 
ically only in the limit of strong nonlinearity can be ob- 
tained in terms of Jacobi elliptic functions for any number 
of particles in the condensate. 

We have also discussed the stability of the stationary 
states under different points of view. The results indicate 



7 



that the soliton-like states, with and without linear coun- 
terpart, are sufficiently stable on the typical time scales 
of a BEC experiment. By voluntarily introducing per- 
turbations of proper intensity, a soliton dynamics could 
also be observed. 
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